import numpy as np
import torch
import matplotlib.pyplot as plt

# 作业4--------------获取正交矩阵
def getMatrix(scopeLeft, scopeRight, n):
    n = int(n)
    # 直接用python自有函数生成
    a = np.random.randint(scopeLeft, scopeRight, n * n).reshape(n, n)
    q, r = np.linalg.qr(a)
    print(q)

# 作业5--------------获取对称正定矩阵symmetrical positive determined matrix
def getSymPosDetMatrix(n):
    # 随机生成矩阵A
    A = torch.tensor(np.random.normal(0, 0.01, (n, n)), dtype=torch.float32)
    B = A.mm(A.t())
    # 确保矩阵对称
    C = B + B.t()
    return C

# 采用梯度下降的方法优化minimizex∈Rn f(x) = 1/2X.TPX + Q.TX + R
# 定义方程
def equation(p, q, r, x):
    return 1/2.*((x.t()).mm(p).mm(x)) + (q.t()).mm(x) + r

# 由公式推导得到的最小的x
def minX(p, q):
    return -1.*(torch.inverse(p).mm(q))

# 定义公式推导出来的梯度
def myGrad(p, q, x):
    return p.mm(x)+q

# 精确线搜索优化学习率
def optimizeLr(a, b, eps, p, q, x, r, lr):
    # 黄金分割点
    t = 0.5 * (5**0.5 - 1)
    # 右搜索点
    lr_r = a + t * (b - a)
    # 左搜索点
    lr_l = a + (1 - t) * (b - a)
    k = 0

    # eps可以接受的搜索区间的值
    while b - a > eps:
        # 如果左点学习率函数值小于右点学习率函数值
        if equation(p, q, r, x-lr_l*myGrad(p, q, x)) < equation(p, q, r, x-lr_r*myGrad(p, q, x)):
            b = lr_r
            a = a
            lr_r = lr_l
            lr_l = a + (1 - t) * (b - a)
        else:
            b = b
            a = lr_l
            lr_l = lr_r
            lr_r = a + t * (b - a)
        k = k + 1
    # 优化好的学习率
    lr = 0.5 * (a + b)
    return lr

# 定义梯度下降+精确直线搜索
# Line Search,这种方法能够保证我们无论从图中哪一点开始梯度下降，都能保证收敛到局部最优值，
# 这种性质也叫做全局收敛性（global convergence）,请一定不要和收敛到全局最优解混淆，
# 全局收敛性只是保证从任意的起点开始都能收敛到局部最优
# 作业6--------------------梯度下降+精确直线搜索”求解上述优化问题
def optimize(n, a, b, eps):
    # 当学习率等于18时，循环5000次可以找到较好的优化解
    # 当学习率小于18时，梯度下降法求最小值时非常缓慢
    # 当学习率为10时，需要循环60000次，才基本找到最小值点
    # 当使用精确线性搜索动态使用学习率时，遍历1000次便可以获得较好的学习结果
    # 注意学习率选取的范围应该取的大一些，不然学习率实际变化并不大
    lr = 2
    err = 1.0
    epochs = 30000
    it = 1
    p = getSymPosDetMatrix(n)
    q = torch.tensor(np.random.normal(0, 0.01, (n, 1)), dtype=torch.float32)
    r = torch.tensor(np.random.normal(0, 0.01, (1, 1)), dtype=torch.float32)
    # 随机变量x,支持梯度下降
    x = torch.tensor(np.random.normal(0, 0.01, (n, 1)), dtype=torch.float32)
    print("p {}".format(p))
    print("q.T {}".format(q.t()))
    print("r {}".format(r))
    print("x.T {}".format(x.t()))

    # 最初的y值
    y_prior = equation(p, q, r, x)
    print("初始值y {}".format(y_prior))
    # 画图像
    curve = [y_prior]
    curve_x1 = [x]
    while err > 1e-2 and it < epochs:
        it += 1
        y = equation(p, q, r, x)
        y.requires_grad_(True)
        # 反向传播
        y.backward()
        # 使用直线搜索法优化学习率
        lr = optimizeLr(a, b, eps, p, q, x, r, lr)
        x = x - myGrad(p, q, x)*lr
        # 清空梯度
        y.grad.zero_()
        # 将x与y放到数组中以绘制图像
        curve.append(y)
        curve_x1.append(x)
    # 公式计算的x，与梯度下降的x
    print("公式x.T {}".format(minX(p, q).t()))
    print("梯度下降x.T {}".format(x.t()))
    print("学习率 {}".format(lr))
    print("实际遍历次数 {}".format(it))
    # ------------y随着遍历次数的增多的变化趋势--------------
    co_x = [i for i in range(it)]
    plt.plot(co_x, curve)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.show()
    # ---------------------{xk的散点图}--------------------
    pic_x = []
    pic_y = []
    for i in range(curve_x1.__len__()):
        pic_x.append(curve_x1[i][0].item())
        pic_y.append(curve_x1[i][1].item())
    plt.scatter(pic_x, pic_y)
    plt.show()

def work7():
    optimize(2, 0, 100000, 0.01)

if __name__ == '__main__':
    # -------------------学习率固定时---------------
    # 作业7
    # optimize(190, 0, 100000000, 0.01)
    optimize(2, 0, 100000, 0.01)
